Skunkware 5/Skunkware 5.iso/lib/calc/cryrand.cal

Text File | 1995-07-17 | 83KB | 2,552 lines

/*
 * cryrand - cryptographically strong pseudo-random number generator library
 */
/*
 * Copyright (c) 1993 by Landon Curt Noll.  All Rights Reserved.
 *
 * Permission to use, copy, modify, and distribute this software and
 * its documentation for any purpose and without fee is hereby granted,
 * provided that the above copyright, this permission notice and text
 * this comment, and the disclaimer below appear in all of the following:
 *
 *    supporting documentation
 *    source copies
 *    source works derived from this source
 *    binaries derived from this source or from derived source
 *
 * LANDON CURT NOLL DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
 * INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS. IN NO
 * EVENT SHALL LANDON CURT NOLL BE LIABLE FOR ANY SPECIAL, INDIRECT OR
 * CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF
 * USE, DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR
 * OTHER TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR
 * PERFORMANCE OF THIS SOFTWARE.
 *
 * chongo was here    /\../\
 */

/*
 * These routines are written in the calc language.  At the time of this
 * notice, calc was at version 1.24.7 (We refer to calc, as in the C
 * program, not the Emacs subsystem).
 *
 * Calc is available by anonymous ftp from ftp.uu.net (137.39.1.9)
 * under the directory /pub/calc.
 *
 * If you can't get calc any other way, Email a request to my Email
 * address as shown below.
 *
 * Comments, suggestions, bug fixes and questions about these routines
 * are welcome.  Send Email to the address given below.
 *
 * Happy bit twiddling,
 *
 *            Landon Curt Noll
 *
 *            [e-mail]
 *            ...!{pyramid,sun,uunet}!hoptoad!chongo
 */

/*
 * AN OVERVIEW OF THE FUNCTIONS:
 *
 * This calc library contains several pseudo-random number generators:
 *
 *   additive 55:
 *
 *    a55rand      - external interface to the additive 55 generator
 *    sa55rand  - seed the additive 55 generator
 *
 *    This is a generator based on the additive 55 generator as
 *    described in Knuth's "The Art of Computer Programming -
 *    Seminumerical Algorithms", vol 2, 2nd edition (1981),
 *    Section 3.2.2, page 27, Algorithm A.
 *
 *    The period and other properties of this generator make it very
 *    useful to 'seed' other generators.
 *
 *    This generator is used by other other generators to produce
 *    various internal values.  In fact, the lower 64 bits of seed
 *    given to other other generators is passed to sa55rand().
 *
 *   shuffle:
 *
 *    shufrand  - generate a pseudo-random value via a shuffle process
 *    sshufrand - seed the shufrand generator
 *    rand      - generate a pseudo-random value over a given range
 *    srand     - seed the random stream generator
 *
 *    This is a generator based on the shuffle generator as described in
 *    Knuth's "The Art of Computer Programming - Seminumerical Algorithms",
 *    vol 2, 2nd edition (1981), Section 3.2.2, page 32, Algorithm B.
 *
 *    The shuffle generator is fast and serves as a fairly good standard
 *    pseudo-random generator.
 *
 *    The shuffle generator is feed random values by the additive 55
 *    generator via a55rand().  Calling a55rand() or sa55rand() will
 *    affect the shuffle generator output.
 *
 *    The rand function is really another interface to the shuffle
 *    generator.  Unlike shufrand(), rand() can return a value of any
 *    given size.  The value returned by rand() may be confined to
 *    either a half open [0,a) (0 <= value < a) or closed interval
 *    [a,b] (a <= value <= b).  By default, a 64 bit value is returned.
 *
 *    Calling srand() simply calls sshufrand() with the same seed.
 *
 *   crypto:
 *
 *    cryrand   - produce a cryptographically strong pseudo-random number
 *    scryrand  - seed the crypto generator
 *    random    - produce a cryptographically strong pseudo-random number
 *            over a given range
 *     srandom   - seed random
 *
 *    This generator is described in the papers:
 *
 *        Blum, Blum, and Shub, "Comparison of Two Pseudorandom Number
 *        Generators", in Chaum, D. et. al., "Advances in Cryptology:
 *        Proceedings Crypto 82", pp. 61-79, Plenum Press, 1983.
 *
 *        "Probabilistic Encryption", Journal of Computer & System
 *        Sciences 28, pp. 270-299.
 *
 *    We also refer to this generator as the 'Blum' generator.
 *
 *    This generator is considered 'strong' in that it passes all
 *    polynomial-time statistical tests.  This means that there
 *    is no statistical test, which runs in polynomial time, that
 *    can distinguish the crypto generator from a truly random source.
 *
 *    The crypto generator is not as fast as most generators, though
 *    it is not painfully slow either.
 *
 *    One may fully seed this generator via scryrand().  Calling
 *    scryrand() with 1 or 3 arguments will result in the additive
 *    55 generator being seeded with the same seed.  Calling
 *     scryrand() with 4 arguments, where the first argument
 *    is >= 0 will also result in the additive 55 generator
 *    being seeded with the same seed.
 *
 *    The random() generator is really another interface to the
 *    crypto generator.  Unlike cryrand(), random() can return a
 *    value confined to either a half open (0 <= value < a) or closed
 *    interval (a <= value <= b).  By default, a 64 bit value is
 *    returned.
 *
 *    Calling srandom() simply calls scryrand(seed).  The additive
 *    55 generator will be seeded with the same seed.
 *
 * All generators come already seeded with precomputed initial constants.
 * Thus, it is not required to seed a generator before using it.
 *
 * Using a seed of '0' will reload generators with their initial states.
 * In the case of scryrand(), lengths of -1 must also be supplied.
 *
 *    sa55rand(0)
 *    sshufrand(0)
 *    srand(0)
 *    scryrand(0,-1,-1)
 *    srandom(0)
 *    scryrand(-1,ip,iq,ir)
 *
 * All of the above 'seed 0' calls are fairly fast.  In fact, passing
 * seeding with a non-zero seed, in the above cases, where seed is
 * not excessively large (many bits long), is also reasonably fast.
 * The 4 arg scryrand() call is fairly fast because no checking is
 * performed on the 'ip', 'iq', or 'ir' values.
 *
 * A call of scryrand(seed,len1,len2), with len1,len2 > 4, (3 args)
 * is a somewhat slow as the length args increase.  This type of
 * call selects 2 primes that are used internally by the crypto
 * generator.  A call of scryrand(seed,ip,iq,ir) where seed >= 0
 * is as slow as the 3 arg case.  See scryrand() for more information.
 *
 * A call of scryrand() (no args) may be used to quickly change the
 * internal state of the crypto and random generators.  Only one major
 * internal crypto generator value (a quadratic residue) is randomly
 * selected via the additive 55 generator.  The other 2 major internal
 * values (the 2 Blum primes) are preserved.  In this form, the additive
 * 55 is not seeded.
 *
 * Calling scryrand() with 1 or 3 args, or calling srandom() will leave
 * the additive 55 generator in a seeded state as if sa55rand(seed) has
 * been called.  Calling scryrand() with 4 args, where the first arg is
 * >0 will also leave the additive 55 generator in the same scryrand(seed)
 * state.  Calling scryrand() with no args will not seed the additive
 * 55 generator before or afterwards, however the additive 55 and shuffle
 * generators will be changed as a side effect of that call.
 *
 * The states of all generators (additive 55, shuffle and crypto) can be
 * saved and restored via the randstate() function.  Saving the state just
 * after * seeding a generator and restoring it later as a very fast way
 * to reseed a generator.
 *
 * As a bonus, the function 'nxtprime' is provided to produce a probable
 * prime number.
 *
 * TRUTH IN ADVERTISING:
 *
 * The word 'probable', in reference to the nxtprime() function, is used
 * because of an extremely extremely small chance that a composite (a
 * non-prime) may be returned.  In no cases will a prime be skipped.
 * For our purposes, this is sufficient as the chance of returning a
 * composite is much smaller than the chance that a hardware glitch
 * will cause nxtprime() to return a bogus result.
 *
 * Another "truth in advertising" issue is the use of the term
 * 'pseudo-random'.  All deterministic generators are pseudo random.
 * This is opposed to true random generators based on some special
 * physical device.
 *
 * The crypto generator is 'pseudo-random'.  There is no statistical
 * test, which runs in polynomial time, that can distinguish the crypto
 * generator from a truly random source.
 *
 * A final "truth in advertising" issue deals with how the magic numbers
 * found in this library were generated.  Detains can be found in the
 * various functions, while a overview can be found in the SOURCE FOR
 * MAGIC NUMBERS section below.
 *
 ****
 *
 * ON THE GENERATORS:
 *
 * The additive 55 generator has a good period, and is fast.  It is
 * reasonable as generators go, though there are better ones available.
 * We use it in seeding the crypto generator as its period and
 * other statistical properties are good enough for our purposes.
 *
 * This shuffle generator has a very good period, and is fast.  It is
 * fairly good as generators go, and is better than the additive 55
 * generator.  Casual direct use of the shuffle generator may be
 * acceptable.  Because of this, the interface to the shuffle generator,
 * but not the additive 55 generator, is advertised when this file is
 * loaded.
 *
 * The shuffle generator functions, shufrand() and rand() use the same
 * seed and tables.  The shuffle generator shuffles the values produced
 * by the additive 55 generator.  Calling or seeding the additive 55
 * generator will affect the output of the shuffle generator.
 *
 * The crypto generator is the best generator in this package.  It
 * produces a cryptographically strong pseudo-random bit sequence.
 * Internally, a fixed number of bits are generated after each
 * generator iteration.  Any unused bits are saved for the next call
 * to the generator.  The crypto generator is not too slow, though
 * seeding the generator from scratch is slow.  Shortcuts and
 * pre-computer seeds have been provided for this reason.  Use of
 * crypto should be more than acceptable for many applications.
 *
 * The crypto seed is in 4 parts, the additive 55 seed (lower 64
 * bits of seed), the shuffle seed (all but the lower 64 bits of
 * seed), and two lengths.  The two lengths specifies the minimum
 * bit size of two primes used internal to the crypto generator.
 * Not specifying the lengths, or using -1 will cause crypto to
 * use the default minimum lengths of 248 and 264 bits, respectively.
 *
 * The random() function just another interface to the crypto
 * generator.  Like rand(), random() provides an interval capability
 * (closed or open) as well as a 64 bit default return value.
 * The random() function as good as crypto, and produces numbers
 * that are equally cryptographically strong.  One may use the
 * seed functions srandom() or scryrand() for either the random()
 * or cryrand() functions.
 *
 * The seed for all of the generators may be of any size.  Only the
 * lower 64 bits of seed affect the additive 55 generator.  Bits
 * beyond the lower 64 bits affect the shuffle generators.  Excessively
 * large values of seed will result in increased memory usage as
 * well as a larger seed time for the shuffle and crypto generators.
 * See REGARDING SEEDS below, for more information.
 *
 * One may save and restore the state of all generators via the
 * randstate() function.
 *
 ****
 *
 * REGARDING SEEDS:
 *
 * Because the generators are interrelated, seeding one generator
 * will directly or indirectly affect the other generators.  Seeding
 * the shuffle generator seeds the additive 55 generator.  Seeding
 * the crypto generator seeds the shuffle generator.
 *
 * The seed of '0' implies that a generator should be seeded back
 * to its initial default state.
 *
 * For the moment, seeds < -1 are reserved for future use.  The
 * value of -1 is used as a special indicator to the fourth form
 * of scryrand(), and it not a real seed.
 *
 * A seed may be of any size.  The additive 55 generator uses only
 * the lower 64 bits, while the shuffle generator uses bytes beyond
 * the lower 64 bits.
 *
 * To help make the generator produced by seed S, significantly
 * different from S+1, seeds are scrambled prior to use.  The
 * function randreseed64() maps [0,2^64) into [0,2^64) in a 1-to-1
 * and onto fashion.
 *
 * The purpose of the randreseed64() is not to add security.  It
 * simply helps remove the human perception of the relationship
 * between the seed and the production of the generator.
 *
 * The randreseed64() process does not reduce the security of the
 * generators.  Every seed is converted into a different unique seed.
 * No seed is ignored or favored.
 *
 * There is no limit on the size of a seed.  On the other hand,
 * extremely large seeds require large tables and long seed times.
 * Using a seed in the range of [2^64, 2^64 * 128!) should be
 * sufficient for most purposes.  An easy way to stay within this
 * range to to use seeds that are between 21 and 215 digits, or 64 to
 * 780 bits long.
 *
 ****
 *
 * SOURCE OF MAGIC NUMBERS:
 *
 * Most of the magic constants used on this library ultimately are
 * based on the Rand book of random numbers.  The Rand book contains
 * 10^6 decimal digits, generated by a physical process.  This book,
 * produced by the Rand corporation in the 1950's is considered
 * a standard against which other generators may be measured.
 *
 * The Rand book of numbers was groups into groups of 20 digits.
 * The first 55 groups < 2^64 were used to initialize add55_init_tbl.
 * The size of 20 digits was used because 2^64 is 20 digits long.
 * The restriction of < 2^64 was used to prevent modulus biasing.
 * (see the note on  modulus biasing in rand()).
 *
 * The additive 55 generator during seeding is used 128 times to help
 * remove the initial seed state from the initial values produced.
 * The loop count of 128 was a power of 2 that permits each of the
 * 55 table entries to be processed at least twice.
 *
 * The function, randreseed64(), uses 4 primes to scramble 64 bits
 * into 64 bits.  These primes were also extracted from the Rand
 * book of numbers.  See sshufrand() for details.
 *
 * The default shuffle table size of 128 entries is the power of 2
 * that is longer than the 100 entries recommended by Knuth for
 * the shuffle algorithm (see the text cited in shufrand()).
 * We use a power of 2 shuffle table length so that the shuffle
 * process can select a table entry from a new additive 55 value
 * by extracting its top most bits.
 *
 * The quadratic residue search performed by cryres() starts at
 * a value that is in the interval [2^sqrpow,n-3), where 2^sqrpow
 * is the first power of 2 > n^(3/4) and n = p*q.  We look at 1009
 * random numbers in this interval for a quadratic residue.  If
 * none are found, sqrpow is decremented by 1 if sqrpow > 2.
 * In practice, decrementing sqrpow never happens.
 *
 * The use of n^(3/4) insures that the quadratic residue is
 * large, but not too large.  We want to avoid residues that are
 * near the square root of 'n', and are nearly 'n' itself (hence
 * the n-3 upper limit) as such residues are trivial or semi-trivial.
 *
 * The '1009' trials is simply an attempt to look for a while before
 * picking a new range.  The number '1009' comes from the first 4
 * digits of the rand book of numbers.
 *
 * Keeping sqrpow > 2 means that we do not look for quadratic
 * residues that are below 4.  This is in keeping with the
 * n-3 in the [2^sqrpow,n-3) interval.
 *
 * The final magic numbers: '248' and '264' are the exponents the
 * largest powers of 2 that are < the two default Blum primes 'p'
 * and 'q' used by the crypto generator.  The values of '248' and
 * '264' implies that the product n=p*q > 2^512.  Each iteration
 * of the crypto generator produces log2(log2(n=p*q)) random bits.
 * The crypto generator is the most efficient when n is slightly >
 * 2^(2^b).  The product n > 2^(2^9)) produces 9 bits as efficiently
 * as possible under the crypto generator process.
 *
 * Not being able to factor 'n=p*q' into 'p' and 'q' does not directly
 * improve the crypto generator.  On the other hand, it can't hurt.
 * The two len values differ slightly to avoid factorization attacks
 * that work on numbers that are a perfect square, or where the two
 * primes are nearly the same.  I elected to have the sizes differ
 * by 3% of the product size.  The difference between '248' and
 * '264', is '16', which is ~3.125% of '512'.  Now 3% of '512' is
 * '~15.36', and the next largest whole number is '16'.
 *
 * The product n=p*q > 2^512 implies a product if at least 155 digits.
 * A product of two primes that is at least 155 digits is slightly
 * beyond Number Theory and computer power of Nov 1992, though this
 * will likely change in the future.
 *
 * Again, the ability (or lack thereof) to factor 'n=p*q' does not
 * directly relate to the strength of the crypto generator.  We
 * selected n=p*q > 2^512 mainly because '512 was a power of 2 and
 * only slightly because it is up in the range where it is difficult
 * to factor.
 *
 ****
 *
 * FOR THE PARANOID:
 *
 * The truly paranoid might suggest that my claims in the MAGIC NUMBERS
 * section are a lie intended to entrap people.  Well they are not, but
 * you need not take my word for it.
 *
 * The random numbers from the Rand book of random numbers can be
 * verified by anyone who obtains the book.  As these numbers were
 * created before I (Landon Curt Noll) was born (you can look up
 * my birth record if you want), I claim to have no possible influence
 * on their generation.
 *
 * There is a very slight chance that the electronic copy of the
 * Rand book that I was given access to differs from the printed text.
 * I am willing to provide access to this electronic copy should
 * anyone wants to compare it to the printed text.
 *
 * One might object to the complexity of the seed scramble/mapping
 * via the randreseed64() function.  The randreseed64() function maps:
 *
 *    1 ==> 4967126403401436567
 *
 * calling srand(1) with the randreseed64() process would be the same
 * as calling srand(4967126403401436567) without it.  No extra security
 * is gained or reduced by using the randreseed64() process.  The meaning
 * of seeds are exchanged, but not lost or favored (used by more than
 * one input seed).
 *
 * One could take issue with my selection of the default sizes '248'
 * and '264'.   As far as I know, 155 digits (512 bits) is beyond the
 * state of the art of Number Theory and Computation as of 01 Sep 92.
 * It will likely be true that 155 digit products of two primes could
 * come within reach in the next few years, but so what?  If you are
 * truly paranoid, why would you want to use the default seed, which
 * is well known?
 *
 * The paranoid today might consider using the lengths of at least '345'
 * and '325' will produce a product of two primes that is 202 digits.
 * (the 2nd and 3rd args of scryrand > 345 & >325 respectively)  Factoring
 * 200+ digit product of two primes is well beyond the current hopes of
 * Number Theory and Computer power, though even this limit may be passed
 * someday.
 *
 * If all the above fails to pacify the truly paranoid, then one may
 * select by some independent means, 2 Blum primes (primes mod 4 == 3,
 * p < q), and a quadratic residue if p*q.  Then by calling:
 *
 *    scryrand(-1, p, q, r)
 *
 * and then calling cryrand() or random(), one may bypass all magic
 * numbers and use the pure generator.
 *
 * Note that randstate() may also be used by the truly paranoid.
 * Even though it holds state for the other generators, their states
 * are independent.
 *
 ****
 *
 * GOALS:
 *
 * The goals of this package are:
 *
 *    all magic numbers are explained
 *
 *        I distrust systems with constants (magic numbers) and tables
 *        that have no justification (e.g., DES).  I believe that I have
 *        done my best to justify all of the magic numbers used.
 *
 *     full documentation
 *
 *        You have this source file, plus background publications,
 *        what more could you ask?
 *
 *    large selection of seeds
 *
 *        Seeds are not limited to a small number of bits.  A seed
 *        may be of any size.
 *
 *    the strength of the generators may be tuned to meet the application need
 *
 *        By using the appropriate seed arguments, one may increase
 *        the strength of the generator to suit the need of the
 *        application.  One does not have just a few levels.
 *
 * This calc lib file is intended for demonstration purposes.  Writing
 * a C program (with possible assembly or libmp assist) would produce
 * a faster generator.
 *
 * Even though I have done my best to implement a good system, you still
 * must use these routines your own risk.
 *
 * Share and enjoy!  :-)
 */


/*
 * These constants are used by all of the generators in various direct and
 * indirect forms.
 */
static cry_seed = 0;            /* master seed */
static cry_mask = 0xffffffffffffffff;    /* 64 bit work mask */


/*
 * Initial magic values for the additive 55 generator - used by sa55rand()
 *
 * These values were generated from the Rand book of random numbers.
 * In groups of 20 digits, we took the first 55 groups < 2^64.
 */
static mat add55_init_tbl[55] = {
    10097325337652013586,    8422689531964509303,
    9376707153831131165,    12807999708015736147,
    12171768336606574717,    2051656926866574818,
    3529647783580834282,    13746700781847540610,
    17468509505804776974,    14385537637435099817,
    14225685144642756788,    11100020401286074697,
    7207317906119690446,    15474452669527079953,
    16868487670307112059,    4493524947524633824,
    13021248927856520106,    15956600001874392423,
    1758753794041921585,    1540354560501451176,
    5335129695612719255,    9973334408846123356,
    2295368703230757546,    15020099946907494138,
    10518216150184876938,    9188200973282539527,
    4220863048338987374,    682273982071453295,
    7706178130835869910,    4618975533122308420,
    397583911260717646,    5686731560708285046,
    10123916228549657560,    1304775865627110086,
    15501295782182641134,    3061180729620744156,
    6958929830512809719,    10850627469959910507,
    13499063195307571839,    6410193623982098952,
    4111084083850807341,    17719042079595449953,
    5462692006544395659,    18288274374963224041,
    8337656769629990836,    7477446061798548911,
    9815931464890815877,    6913451974267278601,
    11883095286301198901,    14974403441045516019,
    14210337129134237821,    12883973436502761184,
    4285013921797415077,    16435915296724552670,
    3742838738308012451
};


/*
 * additive 55 tables - used by a55rand() and sa55rand()
 */
static add55_j = 23;        /* the first walking table pointer */
static add55_k = 54;        /* the second walking table pointer */
static add55_seed64 = 0;    /* lower 64 bits of the reseeded seed */
static mat add55_tbl[55];    /* additive 55 state table */


/*
 * cryobj - cryptographic pseudo-random state object
 */
obj cryobj {                            \
    p,        /* first Blum prime (prime 3 mod 4) */        \
    q,        /* second Blum prime (prime 3 mod 4) */        \
    r,        /* quadratic residue of n=p*q */        \
    exp,    /* used in computing crypto good bits */    \
    left,    /* bits unused from the last cryrand() call */    \
    bitcnt,    /* left contains bitcnt crypto good bits */    \
    a55j,    /* 1st additive 55 table pointer */        \
    a55k,    /* 2nd additive 55 table pointer */        \
    seed,    /* last seed set by sa55rand() or 0 */        \
    shufy,    /* Y (previous a55rand output for shuffle) */    \
    shufsz,    /* size of the shuffle table */            \
    shuftbl,    /* a matrix of shufsz entries */        \
    a55tbl    /* additive 55 generator state table */        \
}


/*
 * initial cryptographic pseudo-random values - used by scryrand()
 *
 * These values are what the crypto generator is initialized with
 * with this library first read.  These values may be reproduced the
 * hard way by calling scryrand(0,248,264) or scryrand(0,-1,-1).
 *
 * We will build up these values a piece at a time to avoid long lines
 * that are difficult to send via Email.
 */
/* p, first Blum prime */
static cryrand_init_p = 0x1aa9e726a735044;
cryrand_init_p <<= 200;
cryrand_init_p |= 0x73b7457c5297122310880fcbfa8d4e38380a00396d533300bb;
/* q, second Blum prime */
static cryrand_init_q = 0xa62ee0481aa12059c3;
cryrand_init_q <<= 200;
cryrand_init_q |= 0x79ef44e72ff58ea70cacabbe9d264a0b16db72117d96f77e17;
/* quadratic residue of n=p*q */
static cryrand_init_r = 0x3d214853f9a5bb4b12f467c9052129a9;
cryrand_init_r <<= 200;
cryrand_init_r |= 0xd215cc6b3c2eae8c7090591b16d8044a886b3c6a58759b1a07;
cryrand_init_r <<= 200;
cryrand_init_r |= 0x2b50a914b42e1b6f9703be86742837c4bc637804c2dc668c5b;

/*
 * cryptographic pseudo-random values - used by cryrand() and scryrand()
 */
/* p, first Blum prime */
static cryrand_p = cryrand_init_p;
/* q, second Blum prime */
static cryrand_q = cryrand_init_q;
/* n = p*q */
static cryrand_n = cryrand_p*cryrand_q;
/* quad residue of n */
static cryrand_r = cryrand_init_r;
/* this cryrand() running exp used in computing crypto good bits */
static cryrand_exp = cryrand_r;
/* good crypto bits unused from the last cryrand() call */
static cryrand_left = 0;
/* the value cryrand_left contains cryrand_bitcnt crypto good bits */
static cryrand_bitcnt = 0;


/*
 * shufrand - shuffle generator constants
 */
static shuf_size = 128;            /* entries in shuffle table */
static shuf_shift = (64-highbit(shuf_size));    /* shift a55 value to get tbl indx */
static shuf_y;                /* Y value (previous output) */
static mat shuf_tbl[shuf_size];        /* shuffle state table */


/*
 * products of consecutive primes - used by nxtprime()
 *
 * We compute these products now to avoid recalculating them on each call.
 * These values help weed out numbers that have a prime factor < 1000.
 */
static nxtprime_pfact10 = pfact(10);
static nxtprime_pfact100 = pfact(100)/nxtprime_pfact10;
static nxtprime_pfact1000 = pfact(1000)/nxtprime_pfact100;


/*
 * a55rand - additive 55 pseudo-random number generator
 *
 * returns:
 *    A number in the half open interval [0,2^64)
 *
 * This function implements the additive 55 pseudo-random number generator.
 *
 * This is a generator based on the additive 55 generator as described in
 * Knuth's "The Art of Computer Programming - Seminumerical Algorithms",
 * vol 2, 2nd edition (1981), Section 3.2.2, page 27, Algorithm A.
 *
 * This function is called from the shuffle generator function shufrand().
 *
 * NOTE: This is NOT Blum's method.  This is used by Blum's method to
 *       generate some internal values.
 *
 * NOTE: If you are looking for a fast generator and do not need a crypto
 *      strong generator, you should not use this generator.  Use of
 *     the shuffle generator functions shufrand() and rand() should
 *     be considered instead.
 */
define
a55rand()
{
    local ret;            /* the pseudo-random number to return */

    /* generate the next value using the additive 55 generator method */
    ret = add55_tbl[add55_k] + add55_tbl[add55_j];
    ret &= cry_mask;
    add55_tbl[add55_k] = ret;

    /* post-process out data with the seed */
    ret = xor(ret, add55_seed64);

    /* step the pointers */
    --add55_j;
    if (add55_j < 0) {
    add55_j = 54;
    }
    --add55_k;
    if (add55_k < 0) {
    add55_k = 54;
    }

    /* return the new pseudo-random number */
    return(ret);
}


/*
 * sa55rand - initialize the additive 55 pseudo-random generator
 *
 * given:
 *    seed
 *
 * returns:
 *    old_seed
 *
 * This function seeds the additive 55 pseudo-random generator.
 *
 * This is a generator based on the additive 55 generator as described in
 * Knuth's "The Art of Computer Programming - Seminumerical Algorithms",
 * vol 2, 2nd edition (1981), Section 3.2.2, page 27, Algorithm A.
 *
 * Unlike Knuth's description, we will let a seed post process our data.
 *
 * We do not apply the seed processing to the additive 55 table
 * data as this would disturb its pseudo-random generator properties.
 * Instead, we xor the output with the low 64 bits of seed.
 *
 * If seed == 0:
 *
 *    This function produces values in exactly the same way
 *    described by Knuth.
 *
 * else:
 *
 *    Each value produced is xor-ed by a 64 bit value.  This value
 *    is the result of randreseed64(rand), which produces a 64
 *    bit value.
 *
 * Anyone comfortable with seed == 0 should also be comfortable with
 * non-zero seeds.  A non-zero seeded sequence will produce a values
 * that have the exact same pseudo-random properties as the algorithm
 * described by Knuth.  I.e., the sequence, while different, is as good
 * (or bad) as the sequence produced by a seed of 0.
 *
 * This function updates both the cry_seed and a55_seed64 global values.
 *
 * This function is called from the shuffle generator seed function sshufrand().
 *
 * NOTE: This is NOT Blum's method.  This is used by Blum's method to
 *       generate some internal values.
 *
 * NOTE: If you are looking for a fast generator and do not need a crypto
 *      strong generator, you should not use this generator.  Use of the
 *     shuffle generator seed functions sshufrand() and srand() should
 *     be considered instead.
 */
define
sa55rand(seed)
{
    local oldseed;        /* previous seed */
    local junk;            /* discards the first few numbers */
    local j;

    /* firewall */
    if (!isint(seed)) {
    quit "bad arg: arg must be an integer";
    }
    if (seed < 0) {
    quit "bad arg: seed < 0 is reserved for future use";
    }

    /* misc table setup */
    oldseed = cry_seed;                /* remember the previous seed */
    cry_seed = seed;                /* save the new seed */
    add55_tbl = add55_init_tbl;    /* reload the table */
    add55_j = 23;               /* reset first walking table pointer */
    add55_k = 54;               /* reset second walking table pointer */

    /* obtain our 64 bit xor seed */
    add55_seed64 = randreseed64(seed);

    /* spin the pseudo-random number generator a while */
    for (j=0; j < 128; ++j) {
    junk = a55rand();
    }

    /* return the old seed */
    return(oldseed);
}


/*
 * shufrand - implement the shuffle pseudo-random number generator
 *
 * returns:
 *    A number in the half open interval [0,2^64)
 *
 * This function implements the shuffle number generator.
 *
 * This is a generator based on the shuffle generator as described in
 * Knuth's "The Art of Computer Programming - Seminumerical Algorithms",
 * vol 2, 2nd edition (1981), Section 3.2.2, page 32, Algorithm B.
 *
 * The function rand() calls this function.
 *
 * NOTE: This is NOT Blum's method.  This is used by Blum's method to
 *       generate some internal values.
 *
 * NOTE: If you do not need a crypto strong pseudo-random generator
 *     this function may very well serve your needs.
 */
define
shufrand()
{
    local j;        /* table index to replace */

    /*
     * obtain a new random value
     * determine the table entry to shuffle
     * shuffle out the value we will return
     */
    shuf_y = shuf_tbl[j = (shuf_y >> shuf_shift)];

    /* shuffle in the new random value */
    shuf_tbl[j] = a55rand();

    /* return the shuffled out value */
    return (shuf_y);
}


/*
 * sshufrand - seed the shuffle pseudo-random generator
 *
 * given:
 *    a seed
 *
 * returns:
 *    the previous seed
 *
 * This function implements the shuffle number generator.
 *
 * This is a generator based on the shuffle generator as described in
 * Knuth's "The Art of Computer Programming - Seminumerical Algorithms",
 * vol 2, 2nd edition (1981), Section 3.2.2, page 32, Algorithm B.
 *
 * The low 64 bits of seed are used by the additive 55 generator.
 * See the sa55rand() function for details.  The remaining bits of seed
 * are used to perform an initial shuffle on the shuffle state table.
 * The size of the seed also determines the size of the shuffle table.
 *
 * The shuffle table size is always a power of 2, and is at least 128
 * entries long.  Let the table size be:
 *
 *    shuf_size = 2^shuf_pow
 *
 * The number of ways one could shuffle that table is:
 *
 *    (2^shuf_pow)!
 *
 * We select the smallest 'shuf_pow' (and thus the size of the shuffle table)
 * such that the following are true:
 *
 *    (2^shuf_pow)! >= (seed / 2^64)    and    2^shuf_pow >= 128
 *
 * Given that we now have the table size of 'shuf_size', we must go about
 * loading the table and shuffling it.
 *
 * Loading is easy, we will generate random values via the additive 55
 * generator (a55rand()) and load them into successive entries.
 *
 * We enumerate the (2^shuf_pow)! values via:
 *
 *    shuf_seed = 2*(U[2] + 3*(U[3] + 4*(U[4] + ...
 *              + (U[shuf_pow-1]*(shuf_pow-1)) ... )))
 *    0 <= U[j] < j
 *
 * We swap the swap table entries shuf_tbl[U[j]] & shuf_tbl[j-1] for all
 * 1 < j < shuf_pow.
 *
 * We will convert 'seed / 2^64' into shuf_seed, by applying the 64 bit
 * scramble function on 64 bit chunks of 'seed / 2^64'.
 *
 * The function srand() calls this function.
 *
 * There is no limit on the size of a seed.  On the other hand,
 * extremely large seeds require large tables and long seed times.
 * Using a seed in the range of [2^64, 2^64 * 128!) should be
 * sufficient for most purposes.  An easy way to stay within this
 * range to to use seeds that are between 21 and 215 digits long, or
 * 64 to 780 bits long.
 *
 * NOTE: This is NOT Blum's method.  This is used by Blum's method to
 *       generate some internal values.
 *
 * NOTE: If you do not need a crypto strong pseudo-random generator
 *     this function may very well serve your needs.
 */
define
sshufrand(seed)
{
    local shuf_pow;        /* power of two - log2(shuf_size) */
    local shuf_seed;        /* seed bits beyond low 64 bits */
    local oldseed;        /* previous seed */
    local shift;        /* shift factor to access 64 bit chunks */
    local swap_indx;        /* what to swap shuf_tbl[0] with */
    local rval;            /* random value form additive 55 generator */
    local j;

    /* firewall */
    if (!isint(seed)) {
    quit "bad arg: seed must be an integer";
    }
    if (seed < 0) {
    quit "bad arg: seed < 0 is reserved for future use";
    }

    /*
     * seed the additive 55 generator
     */
    oldseed = sa55rand(seed);

    /*
     * form the shuffle table size and constants
     */
    shuf_seed = seed >> 64;
    for (shuf_pow = 7; shuf_seed > (j=fact(1<<(shuf_pow))) && shuf_pow < 64; \
     ++shuf_pow) {
    }
    shuf_size = (1 << shuf_pow);
    shuf_shift = 64 - shuf_pow;
    /* reallocate the shuffle table */
    mat shuf_tbl[shuf_size];

    /*
     * scramble the seed above the low 64 bits
     */
    if (shuf_seed > 0) {
    j = 0;
    for (shift=0; shift < highbit(shuf_seed)+1; shift += 64) {
        j |= (randreseed64(shuf_seed >> shift) << shift);
    }
    shuf_seed = j;
    }

    /*
     * load the shuffle table
     */
    for (j=0; j < shuf_size; ++j) {
    shuf_tbl[j] = a55rand();
    }
    shuf_y = a55rand();        /* get the next Y value */

    /*
     * We will shuffle based the process outlined in:
     *
     * Knuth's "The Art of Computer Programming - Seminumerical Algorithms",
     * vol 2, 2nd edition (1981), Section 3.4.2, page 139, Algorithm P.
     *
     * Here, we will let j run over the range [0,shuf_size) instead of
     * [shuf_size,0) as outlined in algorithm P.  We will also generate
     * U values from shuf_seed.
     */
    j = 0;
    while (shuf_seed > 0 && ++j < shuf_size) {

    /* determine what index we will swap with the '0' index */
    quomod(shuf_seed, (j+1), shuf_seed, swap_indx);

    /* swap table entries, if needed */
    if (swap_indx != j) {
        swap(shuf_tbl[j], shuf_tbl[swap_indx]);
    }
    }

    /*
     * run the generator for twice the table size
     */
    for (j=0; j < shuf_size*2; ++j) {
    rval = shufrand();
    }

    /* return the old seed */
    return (oldseed);
}


/*
 * rand - generate a pseudo-random value over a given range via additive 55
 *
 * usage:
 *    rand()         - generate a pseudo-random integer >=0 and < 2^64
 *    rand(a)     - generate a pseudo-random integer >=0 and < a
 *    rand(a,b)    - generate a pseudo-random integer >=a and <= b
 *
 * returns:
 *    a large pseudo-random integer over a give range (see usage)
 *
 * When no arguments are given, a pseudo-random number in the half open
 * interval [0,2^64) is produced.  This form is identical to calling
 * shufrand().
 *
 * When 1 argument is given, a pseudo-random number in the half open interval
 * [0,a) is produced.
 *
 * When 2 arguments are given, a pseudo-random number in the closed interval
 * [a,b] is produced.
 *
 * The input values a and b, if given, must be integers.
 *
 * This generator is simply a different interface to the shuffle generator.
 * calling shufrand(), or seeding via sshufrand() will affect the output
 * of this function.
 *
 * NOTE: Unlike cryrand(), this function does not preserve unused bits for
 *     use by the next function call.
 *
 * NOTE: The Un*x rand() function returns only 16 bit or 31 bits, while we
 *     return a number of any given size (default is 64 bits).
 *
 * NOTE: This is NOT Blum's method.  This is used by Blum's method to
 *       generate some internal values.
 *
 * NOTE: If you do not need a crypto strong pseudo-random generator
 *     this function may very well serve your needs.
 */
define
rand(a,b)
{
    local range;        /* we must generate [0,range) first */
    local offset;        /* what to add to get a adjusted range */
    local ret;            /* pseudo-random value */
    local fullwords;        /* number of full 64 bit chunks in ret */
    local finalmask;        /* mask of bits in final chunk of range */
    local j;

    /*
     * setup and special cases
     */
    /* deal with the rand() case */
    if (isnull(a) && isnull(b)) {
    /* no args means same range as additive 55 generator */
    return(a55rand());
    }
    /* firewall - args, if given must be in integers */
    if (!isint(a) || (!isnull(b) && !isint(b))) {
    quit "bad args: args, if given, must be integers";
    }
    /* convert rand(x) into rand(0,x-1) */
    if (isnull(b)) {
    /* convert call into a closed interval */
    b = a-1;
    a = 0;
    /* firewall - rand(0) should act like rand(0,0) */
    if (b == -1) {
        return(0);
    }
    }
    /* determine the range and offset */
    if (a >= b) {
    /* deal with the case of rand(a,a) */
    if (a == b) {
        /* not very random, but it is true! */
        return(a);
    }
    range = a-b+1;
    offset = b;
    } else {
    /* convert rand(a,b), where a < b into rand(b,a) */
    range = b-a+1;
    offset = a;
    }
    /*
     * At this point, we seek a pseudo-random number [0,range) to which
     * we will add offset to produce a number [offset,range+offset).
     */
    fullwords = highbit(range-1)//64;
    finalmask = (1 << (1+(highbit(range-1)%64)))-1;

    /*
     * loop until we get a value that is in range
     *
     * A note in modulus biasing:
     *
     * We will not fall into the trap of thinking that we can simply take
     * a value mod 'range'.  Consider the case where 'range' is '80'
     * and we are given pseudo-random numbers [0,100).  If we took them
     * mod 80, then the numbers [0,20) would be produced more frequently
     * because the numbers [81,100) mod 80 wrap back into [0,20).
     */
    do {
    /* load up all lower full 64 bit chunks with pseudo-random bits */
    ret = 0;
    for (j=0; j < fullwords; ++j) {
        ret = (ret << 64) + shufrand();
    }

    /* load the highest chunk */
    ret <<= (highbit(finalmask)+1);
    ret |= (shufrand() >> (64-highbit(finalmask)-1));
    } while (ret >= range);

    /*
     * return the adjusted range value
     */
    return(ret+offset);
}


/*
 * srand - seed the pseudo-random additive 55 generator
 *
 * input:
 *    seed
 *
 * returns:
 *    old_seed
 *
 * This function actually seeds the shuffle generator (and indirectly
 * the additive 55 generator used by rand() and a55rand().
 *
 * See sshufrand() and sa55rand() for information about a seed.
 *
 * There is no limit on the size of a seed.  On the other hand,
 * extremely large seeds require large tables and long seed times.
 * Using a seed in the range of [2^64, 2^64 * 128!) should be
 * sufficient for most purposes.  An easy way to stay within this
 * range to to use seeds that are between 21 and 215 digits long, or
 * 64 to 780 bits long.
 *
 * NOTE: This is NOT Blum's method.  This is used by Blum's method to
 *       generate some internal values.
 *
 * NOTE: If you do not need a crypto strong pseudo-random generator
 *     this function may very well serve your needs.
 */
define
srand(seed)
{
    if (!isint(seed)) {
    quit "bad arg: seed must be an integer";
    }
    if (seed < 0) {
    quit "bad arg: seed < 0 is reserved for future use";
    }
    return(sshufrand(seed));
}


/*
 * cryrand - cryptographically strong pseudo-random number generator
 *
 * usage:
 *    cryrand(len)
 *
 * given:
 *    len        number of pseudo-random bits to generate
 *
 * returns:
 *    a cryptographically strong pseudo-random number of len bits
 *
 * Internally, bits are produced log2(log2(n=p*q)) at a time.  If a
 * call to this function does not exhaust all of the collected bits,
 * the unused bits will be saved away and used at a later call.
 * Setting the seed via scryrand() or srandom() will clear out all
 * unused bits.  Thus:
 *
 *    scryrand(0);            <-- restore generator to initial state
 *    cryrand(16);            <-- 16 bits
 *
 * will produce the same value as:
 *
 *    scryrand(0);            <-- restore generator to initial state
 *    cryrand(4)<<12 | cryrand(12);    <-- 4+12 = 16 bits
 *
 * and will produce the same value as:
 *
 *    scryrand(0);            <-- restore generator to initial state
 *    cryrand(3)<<13 | cryrand(7)<<6 | cryrand(6);    <-- 3+7+6 = 16 bits
 *
 * The crypto generator is not as fast as most generators, though it is not
 * painfully slow either.
 *
 * NOTE: This function is the Blum cryptographically strong
 *     pseudo-random number generator!
 */
define
cryrand(len)
{
    local goodbits;    /* the number of good bits generated each pass */
    local goodmask;    /* mask for the low order good bits */
    local randval;    /* pseudo-random value being generated */

    /*
     * firewall
     */
    if (!isint(len) || len < 1) {
    quit "bad arg: len must be an integer > 0";
    }

    /*
     * Determine how many bits may be generated each pass.
     *
     * The result by Alexi et. al., says that the log2(log2(n=p*q))
     * least significant bits are secure, where log2(x) is log base 2.
     */
    goodbits = highbit(highbit(cryrand_n));
    goodmask = (1 << goodbits)-1;

    /*
     * If we have bits left over from the previous call, collect
     * them now.
     */
    if (cryrand_bitcnt > 0) {

    /* case where the left over bits are enough for this call */
    if (len <= cryrand_bitcnt) {

        /* we need only len bits */
        randval = (cryrand_left >> (cryrand_bitcnt-len));

        /* save the unused bits for later use */
        cryrand_left &= ((1 << (cryrand_bitcnt-len))-1);

        /* save away the number of bits that we will not use */
        cryrand_bitcnt -= len;

        /* return our complete result */
        return(randval);

    /* case where we need more than just the left over bits */
    } else {

        /* clear out the number of left over bits */
        len -= cryrand_bitcnt;
        cryrand_bitcnt = 0;

        /* collect all of the left over bits for now */
        randval = cryrand_left;
    }

    /* case where we have no previously left over bits */
    } else {
    randval = 0;
    }

    /*
     * Pump out len cryptographically strong pseudo-random bits,
     * 'goodbits' at a time using Blum's process.
     */
    while (len >= goodbits) {

    /* generate the bits */
    cryrand_exp = (cryrand_exp^2) % cryrand_n;
    randval <<= goodbits;
    randval |= (cryrand_exp & goodmask);

    /* reduce the need count */
    len -= goodbits;
    }

    /* if needed, save the unused bits for later use */
    if (len > 0) {

    /* generate the bits */
    cryrand_exp = (cryrand_exp^2) % cryrand_n;
    randval <<= len;
    randval |= ((cryrand_exp&goodmask) >> (goodbits-len));

    /* save away the number of bits that we will not use */
    cryrand_left = cryrand_exp & ((1 << (goodbits-len))-1);
    cryrand_bitcnt = goodbits-len;
    }

    /*
     * return our pseudo-random bits
     */
     return(randval);
}


/*
 * scryrand - seed the cryptographically strong pseudo-random number generator
 *
 * usage:
 *    scryrand(seed)
 *    scryrand()
 *    scryrand(seed,len1,len2)
 *    scryrand(seed,ip,iq,ir)
 *
 * input:
 *    [seed        pseudo-random seed
 *    [len1 len2]    minimum bit length of the Blum primes 'p' and 'q'
 *            -1 => default lengths
 *    [ip iq ir]    Initial search values for Blum primes 'p', 'q' and
 *            a quadratic residue 'r'
 *
 * returns:
 *    the previous seed
 *
 *
 * This function will seed and setup the generator needed to produce
 * cryptographically strong pseudo-random numbers.  See the function
 * a55rand() and sshufrand() for information about how 'seed' works.
 *
 * The first form of this function are fairly fast if the seed is not
 * excessively large.  The second form is also fairly fast if the internal
 * primes are not too large.  The third form, can take a long time to call.
 * (see below)   The fourth form, if the 'seed' arg is not -1, can take
 * as long as the third form to call.  If the fourth form is called with
 * a 'seed' arg of -1, then it is fairly fast.
 *
 * Calling scryrand() with 1 or 3 args (first and third forms), or
 * calling srandom(), or calling scryrand() with 4 args with the first
 * arg >0, will leave the shuffle generator in a seeded state as if
 * sshufrand(seed) has been called.
 *
 * Calling scryrand() with no args will not seed the shuffle generator,
 * before or afterwards, however the shuffle generator will have been
 * changed as a side effect of that call.
 *
 * Calling scryrand() with 4 args where the first arg is 0 or '-1'
 * will not change the other generators.
 *
 *
 * First form of call:  scryrand(seed)
 *
 * The first form of this function will seed the shuffle generator
 * (via srand).  The default precomputed constants will be used.
 *
 *
 * Second form of call:  scryrand()
 *
 * Only a new quadratic residue of n=p*q is recomputed.  The previous prime
 * values are kept.
 *
 * Unlike the first and second forms of this function, the shuffle
 * generator function is not seeded before or after the call.  The
 * current state is used to generate a new quadratic residue of n=p*q.
 *
 *
 * Third form of call:  scryrand(seed,len1,len2)
 *
 * In the third form, 'len1' and 'len2' guide this function in selecting
 * internally used prime numbers.  The larger the lengths, the longer
 * the time this function will take.  The impact on execution time of
 * cryrand() and random() may also be noticed, but not as much.
 *
 * If a length is '-1', then the default lengths (248 for len1, and 264
 * for len2) are used.  The call scryrand(0,-1,-1) recreates the initial
 * crypto state the slow and hard way.  (use scryrand(0) or srandom(0))
 *
 * This function can take a long time to call given reasonable values
 * of len1 and len2.  On a Pyramid R3000, the time to seed was:
 *
 *    Approx value    digits   seed time
 *      of len1+len2   in n=p*q       in sec
 *    ------------   --------       ------
 *          8           3         0.53
 *         16           5         0.54
 *         32          10         0.79
 *         64          20         1.17
 *        128          39         2.89
 *        200          61         4.68
 *        256          78         7.49
 *        322         100        12.47
 *        464         140        35.56
 *        512         155        53.57
 *        664         200        83.97
 *        830         250       122.93
 *        996         300       242.49
 *       1024         309       295.66
 *       1328         400       663.44
 *       1586         478      2002.10
 *       1660         500      1643.45  (Faster mult/square methods kick in
 *       1992         600      2885.81   in certain cases. Type  help config
 *       2048         617      1578.06   in calc for more details.)
 *
 *     NOTE: The small lengths above are given for comparison
 *           purposes and are NOT recommended for actual use.
 *
 *     NOTE: Generating crypto pseudo-random numbers is MUCH
 *           faster than seeding a crypto generator.
 *
 *     NOTE: This calc lib file is intended for demonstration
 *           purposes.  Writing a C program (with possible assembly
 *           or libmp assist) would produce a faster generator.
 *
 *
 * Fourth form of call:  scryrand(seed,ip,iq,ir)
 *
 * In the fourth form, 'ip', 'iq' and 'ir' serve as initial search
 * values for the two Blum primes 'p' and 'q' and an associated
 * quadratic residue 'r' respectively.  Unlike the 3rd form, where
 * lengths are given, the fourth form allows one to specify minimum
 * search values.
 *
 * The 'seed' value is interpreted as follows:
 *
 *   If seed > 0:
 *
 *    Seed and use the shuffle generator to generate 3 jump values
 *    that are in the range '[0,ip)', '[0,iq)' and '[0,ir)' respectively.
 *    Start searching for legal 'p', 'q' and 'r' values by adding
 *    the jump values to their respective argument values.
 *
 *   If seed == 0:
 *
 *    Start searching for legal 'p', 'q' and 'r' values from
 *    'ip', 'iq' and 'ir' respectively.
 *
 *    This form does not change/seed the other generators.
 *
 *   If seed == -1:
 *
 *    Let 'p' == 'ip', 'q' == 'iq' and 'r' == 'ir'.  Do not check
 *    if the value given are legal Blum primes or an associated
 *    quadratic residue respectively.
 *
 *    This form does not change/seed the other generators.
 *
 *    WARNING: No checks are performed on the args passed.
 *         Passing improper values will likely produce
 *         poor results, or worse!
 *
 *
 * It should be noted that calling scryrand() while using the default
 * primes took only 0.04 seconds.  Calling scryrand(0,-1,-1) took
 * 47.19 seconds.
 *
 * The paranoid, when giving explicit lengths, should keep in mind that
 * len1 and len2 are the largest powers of 2 that are less than the two
 * probable primes ('p' and 'q').  These two primes  will be used
 * internally to cryrand().  For simplicity, we refer to len1 and len2
 * as bit lengths, even though they are actually 1 less then the
 * minimum possible prime length.
 *
 * The actual lengths may exceed the lengths by slightly more than 3%.
 * Furthermore, part of the strength of this generator rests on the
 * difficultly to factor 'p*q'.  Thus one should select 'len1' and 'len2'
 * (from which 'p' and 'q' are selected) such that factoring a 'len1+len2'
 * bit number is difficult.
 *
 * Not being able to factor 'n=p*q' into 'p' and 'q' does not directly
 * improve the crypto generator.  On the other hand, it can't hurt.
 *
 * There is no limit on the size of a seed.  On the other hand,
 * extremely large seeds require large tables and long seed times.
 * Using a seed in the range of [2^64, 2^64 * 128!) should be
 * sufficient for most purposes.  An easy way to stay within this
 * range to to use seeds that are between 21 and 215 digits long, or
 * 64 to 780 bits long.
 *
 * NOTE: This function will clear any internally buffer bits.  See
 *     cryrand() for details.
 *
 * NOTE: This function seeds the Blum cryptographically strong
 *     pseudo-random number generator.
 */
define
scryrand(seed,len1,len2,arg4)
{
    local rval;        /* a temporary pseudo-random value */
    local oldseed;    /* the previous seed */
    local newres;    /* the new quad res */
    local ip;        /* initial Blum prime 'p' search value */
    local iq;        /* initial Blum prime 'q' search value */
    local ir;        /* initial quadratic residue search value */
    local rlim;        /* high quadratic res search value */

    /*
     * firewall - avoid bogus args and very trivial lengths
     */
    /* catch the case of no args - compute a different quadratic residue */
    if (isnull(seed) && isnull(len1) && isnull(len2)) {

    /* generate the next quadratic residue */
    do {
        newres = cryres(cryrand_p, cryrand_q);
    } while (newres == cryrand_r);
    cryrand_r = newres;
    cryrand_exp = cryrand_r;

    /* clear the internal bits */
    cryrand_left = 0;
    cryrand_bitcnt = 0;

    /* return the current seed early */
    return (cry_seed);
    }
    if (!isint(seed)) {
    quit "bad arg: seed arg (1st) must be an integer";
    }
    if (param(0) == 4) {
    if (seed < -1) {
        quit "bad arg: with 4 args: a seed < -1 is reserved for future use";
    }
    } else {
    if (4 && seed < 0) {
        quit "bad arg: a seed arg (1st) < 0 is reserved for future use";
    }
    }

    /*
     * 4 arg case: select or search for 'p', 'q' and 'r' from given values
     */
    if (param(0) == 4) {

    /* set initial values */
    ip = len1;
    iq = len2;
    ir = arg4;

    /*
     * Unless prohibited by a seed if -1, force minimum values on
     * 'ip', 'iq' and 'ir'.
     */
    if (seed >= 0) {
        if (!isint(ip) || ip < 3) {
           /* smallest Blum prime */
           ip = 3;
        }
        if (!isint(iq) || iq < 3) {
           /* smallest Blum prime */
           iq = 3;
        }
        if (!isint(ir) || ir < 4) {
           /* cryres() never selects a value < 4 */
           ir = 4;
        }
    }

    /*
     * Determine our prime search points
     */
    if (seed > 0) {
        /* add in a random value */
        oldseed = srand(seed);
        ip += rand(ip);
        iq += rand(iq);
    }

    /*
     * force ip <= iq
     */
    if (ip > iq) {
        swap(ip, iq);
    }

    /*
     * find the first Blum prime 'p'
     */
    if (seed >= 0) {
        cryrand_p = nxtprime(ip,3,4);
    } else {
        /* seed == -1: assume 'ip' is a Blum prime */
        cryrand_p = ip;
    }

    /*
     * find the 2nd Blum prime 'q' > 'p', if needed
     */
    if (seed >= 0) {
        if (iq <= cryrand_p) {
        iq = cryrand_p + 2;
        }
        cryrand_q = nxtprime(iq,3,4);
    } else {
        /* seed == -1: assume 'iq' is a Blum prime */
        cryrand_q = iq;
    }

    /* remember our p*q Blum prime product as well */
    cryrand_n = cryrand_p*cryrand_q;

    /*
     * search for a quadratic residue
     */
    if (seed >= 0) {

        /* add in a random value to 'ir' if seeded */
        if (seed > 0) {
        ir += rand(ir);
        }

        /*
         * increment until we find a quadratic value
         *
         * p is prime and J(r,p) == 1  ==>  r is a quadratic residue of p
         * q is prime and J(q,p) == 1  ==>  r is a quadratic residue of q
         *
         * J(r,p)*J(r,q) == J(r,p*q)
         * J(r,p) and J(q,p) == 1       ==>  J(r,p*q) == 1
         * J(r,p*q) == 1           ==>  r is a quadratic residue of p*q
         *
         * Try to avoid trivial quadratic residues by forcing the search
         * over the interval [4, n-3).
         */
        rlim = cryrand_n-3;
        /* if ir is too big, drop down to the bottom of the range */
        if (ir >= rlim) {
        ir = 4;
        }
        while (jacobi(ir,cryrand_p) != 1 || jacobi(ir,cryrand_q) != 1) {
        /* if ir is too big, drop down to the bottom of the range */
        if (++ir >= rlim) {
            ir = 4;
        }
        }
    }
    cryrand_r = ir;

    /*
     * clear the previously unused cryrand bits & other misc setup
     */
    cryrand_left = 0;
    cryrand_bitcnt = 0;
    cryrand_exp = cryrand_r;

    /*
     * reseed the generator, if needed
     *
     * The crypto generator no longer needs the additive 55 and shuffle
     * generators.  We however, restore the additive 55 and shuffle
     * generators back to its seeded state in order to be sure that it
     * will be left in the same state.
     *
     * This will make more reproducible, calls to the additive 55 and
     * shuffle generator; or more reproducible, calls to this function
     * without args.
     */
    if (seed > 0) {
        ir = srand(seed);   /* ignore this return value */
        return(oldseed);
    } else {
        /* no seed change, return the current seed */
        return (cry_seed);
    }
    }

    /*
     * If not the 4 arg case:
     *
     * convert explicit -1 args into default values
     * convert missing args into -1 values (use precomputed tables)
     */
    if ((isint(len1) && len1 != -1 && len1 < 4) ||
    (isint(len2) && len2 != -1 && len2 < 4) ||
    (!isint(len1) && isint(len2)) ||
    (isint(len1) && !isint(len2))) {
    quit "bad args: 2 & 3: if 2nd is given, must be -1 or ints > 3";
    }
    if (isint(len1) && len1 == -1) {
    len1 = 248;    /* default len1 value */
    }
    if (isint(len2) && len2 == -1) {
    len2 = 264;    /* default len2 value */
    }
    if (!isint(len1) && !isint(len2)) {
    /* from here down, -1 means use precomputed values */
    len1 = -1;
    len2 = -1;
    }

    /*
     * force len1 <= len2
     */
    if (len1 > len2) {
    swap(len1, len2);
    }

    /*
     * seed the generator
     */
    oldseed = srand(seed);

    /*
     * generate p and q Blum primes
     *
     * The Blum process requires the primes to be of the form 3 mod 4.
     * We also generate n=p*q for future reference.
     *
     * We make sure that the lengths are the minimum lengths possible.
     * We want some range to select a random prime from, so we
     * go at least 3 bits higher, and as much as 3% plus 3 bits
     * higher.  Since the section is a random, how high really
     * does not matter that much, but we want to avoid going to
     * an extreme to keep the execution time from getting too long.
     *
     * Finally, we generate a quadratic residue of n=p*q.
     */
    if (len1 > 0) {
    /* generate a pseudo-random prime ~len1 bits long */
    rval = rand(2^(len1-1), 2^((int(len1*1.03))+3));
    cryrand_p = nxtprime(rval,3,4);
    } else {
    /* use precomputed 'p' value */
    cryrand_p = cryrand_init_p;
    }
    if (len2 > 0) {
    /* generate a pseudo-random prime ~len1 bits long */
    rval = rand(2^(len2-1), 2^((int(len2*1.03))+3));
    cryrand_q = nxtprime(rval,3,4);
    } else {
    /* use precomputed 'q' value */
    cryrand_q = cryrand_init_q;
    }

    /*
     * find the quadratic residue
     */
    cryrand_n = cryrand_p*cryrand_q;
    cryrand_r = cryres(cryrand_p, cryrand_q);

    /*
     * clear the previously unused cryrand bits & other misc setup
     */
    cryrand_left = 0;
    cryrand_bitcnt = 0;
    cryrand_exp = cryrand_r;

    /*
     * reseed the generator
     *
     * The crypto generator no longer needs the additive 55 and shuffle
     * generators.  We however, restore the additive 55 and shuffle
     * generators back to its seeded state in order to be sure that it
     * will be left in the same state.
     *
     * This will make more reproducible, calls to the additive 55 and
     * shuffle generator; or more reproducible, calls to this function
     * without args.
     */
    /* we do not care about this old seed */
    rval = srand(seed);

    /*
     * return the old seed
     */
    return(oldseed);
}


/*
 * random - a cryptographically strong pseudo-random number generator
 *
 * usage:
 *    random()     - generate a pseudo-random integer >=0 and < 2^64
 *    random(a)     - generate a pseudo-random integer >=0 and < a
 *    random(a,b)    - generate a pseudo-random integer >=a and <= b
 *
 * returns:
 *    a large cryptographically strong pseudo-random number  (see usage)
 *
 * This function is just another interface to the crypto generator.
 * (see the cryrand() function).
 *
 * When no arguments are given, a pseudo-random number in the half open
 * interval [0,2^64) is produced.  This form is identical to calling
 * cryrand(64).
 *
 * When 1 argument is given, a pseudo-random number in the half open interval
 * [0,a) is produced.
 *
 * When 2 arguments are given, a pseudo-random number in the closed interval
 * [a,b] is produced.
 *
 * This generator uses the crypto to return a large pseudo-random number.
 *
 * The input values a and b, if given, must be integers.
 *
 * Internally, bits are produced log2(log2(n=p*q)) at a time.  If a
 * call to this function does not exhaust all of the collected bits,
 * the unused bits will be saved away and used at a later call.
 * Setting the seed via scryrand(), srandom() or cryrand(len,1)
 * will clear out all unused bits.
 *
 * NOTE: The BSD random() function returns only 31 bits, while we return 64.
 *
 * NOTE: This function is the Blum cryptographically strong
 *     pseudo-random number generator!
 */
define
random(a,b)
{
    local range;        /* we must generate [0,range) first */
    local offset;        /* what to add to get a adjusted range */
    local rangebits;        /* the number of bits in range */
    local ret;            /* pseudo-random bit value */

    /*
     * setup and special cases
     */
    /* deal with the rand() case */
    if (isnull(a) && isnull(b)) {
    /* no args means return 64 bits */
    return(cryrand(64));
    }
    /* firewall - args, if given must be in integers */
    if (!isint(a) || (!isnull(b) && !isint(b))) {
    quit "bad args: args, if given, must be integers";
    }
    /* convert rand(x) into rand(0,x-1) */
    if (isnull(b)) {
    /* convert call into a closed interval */
    b = a-1;
    a = 0;
    /* firewall - rand(0) should act like rand(0,0) */
    if (b == -1) {
        return(0);
    }
    }
    /* determine the range and offset */
    if (a >= b) {
    /* deal with the case of rand(a,a) */
    if (a == b) {
        /* not very random, but it is true! */
        return(a);
    }
    range = a-b+1;
    offset = b;
    } else {
    /* convert random(a,b), where a<b, into random(b,a) */
    range = b-a+1;
    offset = a;
    }
    rangebits = highbit(range-1)+1;
    /*
     * At this point, we seek a pseudo-random number [0,range) to which
     * we will add offset to produce a number [offset,range+offset).
     */

    /*
     * loop until we get a value that is in range
     *
     * We will obtain pseudo-random values over the range [0,2^rangebits)
     * where 2^rangebits >= range and 2^(rangebits-1) < range.  We
     * will ignore any results that are > the range that we want.
     *
     * A note in modulus biasing:
     *
     * We will not fall into the trap of thinking that we can simply take
     * a value mod 'range'.  Consider the case where 'range' is '80'
     * and we are given pseudo-random numbers [0,100).  If we took them
     * mod 80, then the numbers [0,20) would be produced more often
     * because the numbers [81,100) mod 80 wrap back into [0,20).
     */
    do {
    /* obtain a pseudo-random value */
    ret = cryrand(rangebits);
    } while (ret >= range);

    /*
     * return the adjusted range value
     */
    return(ret+offset);
}


/*
 * srandom - seed the cryptographically strong pseudo-random number generator
 *
 * given:
 *    seed    a random number seed
 *
 * returns:
 *      the previous seed
 *
 * This function is just another interface to the crypto generator.
 * (see the scryrand() function).
 *
 * This function makes indirect use of the additive 55 and shuffle
 * generator.
 *
 * There is no limit on the size of a seed.  On the other hand,
 * extremely large seeds require large tables and long seed times.
 * Using a seed in the range of [2^64, 2^64 * 128!) should be
 * sufficient for most purposes.  An easy way to stay within this
 * range to to use seeds that are between 21 and 215 digits long, or
 * 64 to 780 bits long.
 *
 * NOTE: Calling this function will clear any internally buffer bits.
 *     See cryrand() for details.
 *
 * NOTE: This function seeds the Blum cryptographically strong
 *     pseudo-random number generator!
 */
define
srandom(seed)
{
    if (!isint(seed)) {
    quit "bad arg: seed must be an integer";
    }
    if (seed < 0) {
    quit "bad arg: seed < 0 is reserved for future use";
    }
    return(scryrand(seed));
}


/*
 * randstate - set/get the state of all of the generators
 *
 * usage:
 *    randstate()    return the current state
 *    randstate(0)    return the previous state, set the default state
 *    randstate(cobj)    return the previous state, set a new state
 *
 * In the first form: randstate()
 *
 *    This function returns an cryobj object containing information
 *    about the current state of all of the generators.
 *
 * In the second form: randstate(0)
 *
 *    This function sets all of the generators to the default initial
 *    state (i.e., the state when this library was loaded).
 *
 *    This function returns an cryobj object containing information
 *    about the previous state of all of the generators.
 *
 * In the third form: randstate(cobj)
 *
 *    This function sets all of the generators to the state as found
 *    in the cryobj object.
 *
 *    This function returns an cryobj object containing information
 *    about the previous state of all of the generators.
 *
 * This function may be used to save and restore cryrand() & random()
 * generator states.  For example:
 *
 *    state = randstate()        <-- save the current state
 *    random()            <-- print the next 64 bits
 *    randstate(state)        <-- restore previous state
 *    random()            <-- print the same 64 bits
 *
 * One may quickly reseed a generator.  For example:
 *
 *    srandom(1,330,350)        <-- seed the generator
 *    seed1state = randstate()    <-- remember this 1st seeded state
 *    random()            <-- print 1st 64 bits seed 1 generator
 *    srandom(2,331,351)        <-- seed the generator again
 *    seed2state = randstate()    <-- remember this 2nd seeded state
 *    random()            <-- print 1st 64 bits seed 2 generator
 *
 *    randstate(seed1state)        <-- reseed to the 1st seeded state
 *    random()            <-- reprint 1st 64 bits seed 1 generator
 *    randstate(seed2state)        <-- reseed to the 2nd seeded state
 *    random()            <-- reprint 1st 64 bits seed 2 generator
 *
 *    oldstate = randstate(0)        <-- seed to the default generator
 *    random()            <-- print 1st 64 bits from default
 *    a55rand()            <-- print 1st 64 bits a55 generator
 *    prevstate = randstate(oldstate)    <-- restore seed 2 generator
 *    random()            <-- print 2nd 64 bits seed 2 generator
 *    randstate(prevstate)        <-- restore def generator in progress
 *    random()            <-- print 2nd 64 bits default generator
 *    a55rand()            <-- print 2nd 64 bits a55 generator
 *
 * given:
 *    cobj    if a cryobj object, use that object to set the current state
 *        if 0, set to the default state
 *
 * return:
 *    return the state of the crypto pseudo-random number generator in
 *           the form of an cryobj object, as it was prior to this call
 *
 * NOTE: No checking is performed on the data the 3rd form (cryobj object
 *     arg) is used.  The user must ensure that the arg represents a valid
 *     generator state.
 *
 * NOTE: When using the second form (passing an integer arg), only 0 is
 *     defined.  All other integer values are reserved for future use.
 */
define
randstate(arg)
{
    /* declare our objects */
    local obj cryobj x;        /* firewall comparator */
    local obj cryobj prev;    /* previous states of the generators */
    local junk;            /* dummy holder of random values */

    /* firewall */
    if (!isint(arg) && !istype(arg,x) && !isnull(arg)) {
    quit "bad arg: argument must be integer, an cryobj object or missing";
    }
    if (isint(arg) && arg != 0) {
        quit "bad arg:  non-zero integer arguments are reserved for future use";
    }

    /*
     * save the current state
     */
    prev.p = cryrand_p;
    prev.q = cryrand_q;
    prev.r = cryrand_r;
    prev.exp = cryrand_exp;
    prev.left = cryrand_left;
    prev.bitcnt = cryrand_bitcnt;
    prev.a55j = add55_j;
    prev.a55k = add55_k;
    prev.seed = cry_seed;
    prev.shufy = shuf_y;
    prev.shufsz = shuf_size;
    prev.shuftbl = shuf_tbl;
    prev.a55tbl = add55_tbl;
    if (isnull(x)) {
    /* if no args, just return current state */
    return (prev);
    }

    /*
     * deal with the cryobj arg - set the state
     */
    if (istype(arg, x)) {
    /* set the state from this object */
    cryrand_p = arg.p;
    cryrand_q = arg.q;
    cryrand_n = cryrand_p*cryrand_q;
    cryrand_r = arg.r;
    cryrand_exp = arg.exp;
    cryrand_left = arg.left;
    cryrand_bitcnt = arg.bitcnt;
    add55_j = arg.a55j;
    add55_k = arg.a55k;
    cry_seed = arg.seed;
    add55_seed64 = randreseed64(cry_seed);
    shuf_y = arg.shufy;
    shuf_size = arg.shufsz;
    shuf_shift = (64-highbit(shuf_size));
    mat shuf_tbl[shuf_size];
    shuf_tbl = arg.shuftbl;
    add55_tbl = arg.a55tbl;

    /*
     * deal with the 0 integer arg - set the default initial state
     */
    } else if (isint(arg) && arg == 0) {
    cryrand_p = cryrand_init_p;
    cryrand_q = cryrand_init_q;
    cryrand_n = cryrand_p * cryrand_q;
    cryrand_r = cryrand_init_r;
    cryrand_exp = cryrand_r;
    cryrand_left = 0;
    cryrand_bitcnt = 0;
    cry_seed = 0;
    cry_seed = sshufrand(0);
    }

    /*
     * return the previous state
     */
    return (prev);
}


/*
 * nxtprime - find a probable prime >= n_arg
 *
 * usage:
 *    nxtprime(n_arg)
 *    nxtprime(n_arg, modval, modulus)
 *
 * given:
 *    n_arg            lower bound of the search
 *    [modval modulus]    if given, look for numbers mod modulus == modval
 *
 * returns:
 *    A number is that is very likely prime.
 *
 * In the first case 'nxtprime(n_arg)', this function returns a probable
 * prime >= n_arg.  In the second case 'nxtprime(n_arg, v, u)', this
 * function returns a probable prime >= n_arg that is also == v mod u.
 *
 * This function will not skip over a prime, through there is a
 * extremely unlikely chance that it will return a composite.
 * The odds that a number returned by this function is not prime
 * are 1 in 4^50.  The failure rate of this function is many orders
 * or magnitude lower than the failure rate due to a hardware error.
 *
 * NOTE: This function can take a long time, given a large value of n_arg.
 */
define
nxtprime(n_arg, modval, modulus)
{
    local modgcd;        /* gcd(modulus,modval) */
    local n;            /* value >= n_arg that is being tested */
    local j;

    /*
     * firewall
     */
    if (!isint(n_arg)) {
    quit "bad args: 1st arg must be an integer";
    }
    if (!isnull(modulus) && !isint(modval)) {
    quit "bad args: 3rd arg, if 2nd arg is given, must be an integer";
    }
    if (!isnull(modulus) && (!isint(modulus) || modulus <= 0)) {
    quit "bad args: 3nd arg, if given, must be an integer > 0";
    }

    /*
     * get values < 3 out of the way
     */
    n = n_arg;
    if (n < 3) {
    /* get the even prime out of the way, if possible */
    if (isnull(modulus) ||
        modulus == 1 ||
        (n%modulus == modval%modulus)) {
        /*
         * 2 is the greatest odd prime, because
         * 2 is the least even prime  :-)
         */
        return(2);
    }
    /* we have eliminated everything < 3 */
    n = 3;
    }

    /*
     * convert nxtprime(n) to nxtprime(n,1,2)
     * convert nxtprime(n,x,1) to nxtprime(n,1,2)
     * convert nxtprime(n,a,b) to nxtprime(n,a mod b,b)
     */
    if (isnull(modulus) || modulus < 2) {
    modulus = 2;
    modval = 1;
    }
    modval %= modulus;

    /*
     * catch cases where no more primes == 'modval' mod 'modulus' exist
     */
    modgcd = gcd(modval,modulus);
    if (modgcd > 1) {

    /* if beyond the modgcd, then no primes can exist */
        if (n > modgcd) {
        print "n_arg:",n_arg,"  modval:",modval,"  modulus:",modulus;
        quit "no such prime of that form exists";
    }

    /* else n <= modgcd, then our only chance is if modgcd is prime */
    /* reject if modgcd has an obvious prime factor */
    if (modgcd > 10 && gcd(modgcd,nxtprime_pfact10) != 1) {
        print "n_arg:",n_arg,"  modval:",modval,"  modulus:",modulus;
        quit "no such prime of that form exists";
    }
    if (modgcd > 100 && gcd(modgcd,nxtprime_pfact100) != 1) {
        print "n_arg:",n_arg,"  modval:",modval,"  modulus:",modulus;
        quit "no such prime of that form exists";
    }
    if (modgcd > 1000 && gcd(modgcd,nxtprime_pfact1000) != 1) {
        print "n_arg:",n_arg,"  modval:",modval,"  modulus:",modulus;
        quit "no such prime of that form exists";
    }

    /* do 50 probable prime tests, for a 1 in 4^50 false prime rate */
    if (!ptest(modgcd,50)) {
        print "n_arg:",n_arg,"  modval:",modval,"  modulus:",modulus;
        quit "no such prime of that form exists";
    }

    /* modgcd is the only prime >= n */
    return(modgcd);
    }

    /*
     * bump n up to the next possible case
     *
     * n will be an odd number == 'modval' mod 'modulus'
     */
    if (n%modulus != modval) {
    j = n - (n%modulus) + modval;
    if (j < n) {
        n = j+modulus;
    } else {
        n = j;
    }
    }
    if (n%2 == 0) {
    n += modulus;
    }

    /* look for a prime */
    n = n-modulus;
    do {
    /* try the next integer */
    n = n+modulus;

    /* reject if it has an obvious prime factor */
    if (n > 10 && gcd(n,nxtprime_pfact10) != 1) {
        continue;
    }
    if (n > 100 && gcd(n,nxtprime_pfact100) != 1) {
        continue;
    }
    if (n > 1000 && gcd(n,nxtprime_pfact1000) != 1) {
        continue;
    }

    /* do 50 probable prime tests */
    if (!ptest(n,50)) {
        continue;
    }

    /* n is very likely a prime number */
    return(n);

    } while (1);
}


/*
 * cryobj - how to initialize a cryobj object
 *
 * given:
 *    p        first Blum prime (prime 3 mod 4)
 *    q        second Blum prime (prime 3 mod 4)
 *    r        quadratic residue of n=p*q
 *    exp        used in computing crypto good bits
 *    left        bits unused from the last cryrand() call
 *    bitcnt        left contains bitcnt crypto good bits
 *    a55j        1st additive 55 table pointer
 *    a55k        2nd additive 55 table pointer
 *    seed        last seed set by sa55rand() or 0
 *    shufy        Y (previous a55rand() output for shuffle)
 *    shufsz        size of the shuffle table
 *    shuftbl        a matrix of shufsz entries
 *    a55tbl        additive 55 generator state table
 *
 * return:
 *    an cryobj object
 *
 * NOTE: This function, by convention, returns an cryobj object.
 */
define
cryobj(p,q,r,exp,left,bitcnt,a55j,a55k,seed,shufy,shufsz,shuftbl,a55tbl)
{
    /* declare our objects */
    local obj cryobj x;

    /* firewall */
    if (!isint(p) || !isint(q) || !isint(r) || !isint(exp) || \
    !isint(left) || !isint(bitcnt) || !isint(a55j) || \
    !isint(a55k) || !isint(seed) || !isint(shufy) || !isint(shufsz)) {
    quit "bad args: first 11 args must be integers";
    }
    if (!ismat(shuftbl) || matdim(shuftbl) != 1 || \
    matmin(shuftbl,1) != 0 || matmax(shuftbl,1) != shuf_size-1) {
    quit "bad arg: 12th is not a mat[0:shuf_size-1]";
    }
    if (!ismat(a55tbl) || matdim(a55tbl) != 1 || \
    matmin(a55tbl,1) != 0 || matmax(a55tbl,1) != 54) {
    quit "bad arg: 13th is not a mat[0:54]";
    }

    /* initialize object with default startup values */
    x.p = p;
    x.q = q;
    x.r = r;
    x.exp = exp;
    x.left = left;
    x.bitcnt = bitcnt;
    x.a55j = a55j;
    x.a55k = a55k;
    x.seed = seed;
    x.shufy = shuf_y;
    x.shufsz = shuf_size;
    x.shuftbl = shuf_tbl;
    x.a55tbl = a55tbl;

    /* return the initialized object */
    return (x);
}


/*
 * cryobj_print - print the value of a cryobj object
 *
 * usage:
 *    a    an cryobj object
 *
 * NOTE: This function is called automatically when an cryobj object
 *       is displayed.
 */
define
cryobj_print(a)
{
    /* declare our objects */
    local obj cryobj x;        /* firewall comparator */

    /* firewall */
    if (!istype(a, x)) {
    quit "bad arg: arg is not an cryobj object";
    }
    if (!ismat(a.shuftbl) || matdim(a.shuftbl) != 1 || \
    matmin(a.shuftbl,1) != 0 || matmax(a.shuftbl,1) != a.shufsz-1) {
    quit "bad arg: 12th is not a mat[0:shuf_size-1]";
    }
    if (!ismat(a.a55tbl) || matdim(a.a55tbl) != 1 || \
    matmin(a.a55tbl,1) != 0 || matmax(a.a55tbl,1) != 54) {
    quit "bad arg: 13th is not a mat[0:54]";
    }

    /* print the value */
    print "cryobj(" : a.p : "," : a.q : "," : a.r : "," : a.exp : "," : \
      a.left : "," : a.bitcnt : "," : a.a55j : "," : a.a55k : "," : \
      a.seed : "," : a.shufy : "," : a.shufsz : \
      ",[" : a.shuftbl[0] : "," : a.shuftbl[1] : "," : \
      a.shuftbl[2] : ",...," : a.shuftbl[52] : "," : \
      a.shuftbl[53] : "," : a.shuftbl[54] : "]" : \
      ",[" : a.a55tbl[0] : "," : a.a55tbl[1] : "," : \
      a.a55tbl[2] : ",...," : a.a55tbl[52] : "," : \
      a.a55tbl[53] : "," : a.a55tbl[54] : "])" : ;
}


/*
 * cryres - find a pseudo-random quadratic residue for scryrand() and cryrand()
 *
 * given:
 *    p    first prime
 *    q    second prime
 *
 * returns:
 *    a number that is a quadratic residue of n=p*q
 *
 * This function is returns the pseudo-random quadratic residue of
 * the product of two primes.  Normally this function is called
 * only by the crypto generator.
 *
 * NOTE: No check is made to ensure that the values passed are primes.
 *
 * NOTE: This is NOT Blum's method.  This is used by Blum's method to
 *       generate some internal values.
 */
define
cryres(p,q)
{
    local rval;        /* a temporary pseudo-random value */
    local sqrpow;    /* a power of 2 starting > n^(3/4) */
    local quadres;    /* 0 or a quadratic residue of n = p*q */
    local n;        /* n=p*q */
    local j;

    /*
     * firewall
     */
    if (!isint(p) || !isint(q) || p < 3 || q < 3) {
    quit "bad args: p and q must be integers > 2";
    }

    /*
     * find a pseudo-random quadratic residue of n = p*q
     *
     * To find a pseudo-random quadratic residue, we will generate
     * pseudo-random numbers to try in the interval [2^sqrpow,n-3),
     * where 2^sqrpow is the first power of 2 > n^(3/4).  This range
     * helps us avoid selecting trivial residues abs(quadres mod n=p*q)
     * is tiny.  We will never select a residue < 4.
     *
     * When we fail to find a quadratic residue after 1009 tries, we will
     * drop sqrpow by 1 and start at another pseudo-random location.
     *
     * It is very unlikely that we will need to search more than a
     * few numbers to find a quadratic residue of n = p*q.
     *
     * We will reject any quadratic residue that is a square of
     * itself, mod n=p*q.
     */
    n = p*q;
    quadres = 0;
    sqrpow = highbit(n)//4*3;
    do {
    /* generate a pseudo-random starting point */
    rval = rand(2^sqrpow, n-4);

    /* look for 1009 times or until >= cryrand_n */
    for (j=0; j < 1009; ++j, ++rval) {

        /*
         * check if we have a quadratic residue of n = p*q
         *
         * p is prime and J(r,p) == 1  ==>  r is a quadratic residue of p
         * q is prime and J(q,p) == 1  ==>  r is a quadratic residue of q
         *
         * J(r,p)*J(r,q) == J(r,p*q)
         * J(r,p) and J(q,p) == 1       ==>  J(r,p*q) == 1
         * J(r,p*q) == 1           ==>  r is a quadratic residue of p*q
         */
        if (jacobi(rval,p) == 1 && jacobi(rval,q) == 1) {
        quadres = rval;
        break;
        }
    }

    /* setup for a new pseudo-random starting point if we found nothing */
    if (quadres == 0) {
        /* reduce sqrpow if reasonable */
        if (sqrpow > 2) {
        --sqrpow;
        }
    }
    } while (quadres == 0 || ((quadres^2)%n) == quadres);

    /*
     * return the quadratic residue of n = p*q
     */
    return (quadres);
}


/*
 * randreseed64 - scramble a 64 bit seed
 *
 * given:
 *    a 64 bit seed
 *
 * returns:
 *    a 64 scrambled seed, or 0 if seed was 0
 *
 * It is 'nice' when a seed of "n" produces a 'significantly different'
 * sequence than a seed of "n+1".  Generators, by convention, assign
 * special significance to the seed of '0'.  It is an unfortunate that
 * people often pick small seed values, particularly when large seed
 * are of significance to the generators found in this file.
 *
 * If 'seed' is 0, then 0 is returned.  If 'seed' is non-zero, we will
 * produce a different and unique new scrambled 'seed'.  This scrambling
 * will effectively eliminate the human factors and perceptions that
 * are noted above.
 *
 * It should be noted that the purpose of this process to scramble a seed
 * ONLY.  We do not care if these generators produce good random numbers.
 * We only want to help eliminate the human factors and perceptions
 * noted above.
 *
 * This function scrambles the low 64 bits of a seed, by mapping [0,2^64)
 * into [0,2^64).  This map is one-to-one and onto.  Mapping is performed
 * using  a linear congruence generator of the form:
 *
 *        X1 <-- (a*X0 + c) mod m
 *
 * The generator are based on the linear congruential generators found in
 * Knuth's "The Art of Computer Programming - Seminumerical Algorithms",
 * vol 2, 2nd edition (1981), Section 3.6, pages 170-171.
 *
 * Because we process 64 bits we will take:
 *
 *        m = 2^64            (based on note ii)
 *
 * We will scan the Rand book of numbers to look for an 'a' such that:
 *
 *        a mod 8 == 5            (based on note iii)
 *        0.01*m < a < 0.99*m        (based on note iv)
 *        0.01*2^64 < a < 0.99*2^64
 *
 * To help keep the generators independent, we want:
 *
 *        a is prime
 *
 * The choice of an adder 'c' is considered immaterial according (based
 * in note v).  Knuth suggests 'c==1' or 'c==a'.  We elect to select 'c'
 * using the same process as we used to select 'a'.  The choice is
 * 'immaterial' after all, and as long as:
 *
 *        gcd(c, m) == 1        (based on note v)
 *        gcd(c, 2^64) == 1
 *
 * the concerns are met.   It can be shown that if we have:
 *
 *        gcd(a, c) == 1
 *
 * then the adders and multipliers will be more independent.
 *
 * We will obtain the values 'a' and 'c for our generator from the
 * Rand book of numbers.  Because m=2^64 is 20 decimal digits long, we
 * will search the Rand book of numbers 20 at a time.  We will skip any
 * of the 55 values that were used to initialize the additive 55 generators.
 * The values obtained from the Rand book are:
 *
 *        a = 6316878969928993981
 *        c = 1363042948800878693
 *
 * As we stated before, we must map 0 ==> 0.  The linear congruence
 * generator would normally map as follows:
 *
 *    0 ==> 1363042948800878693    (0 ==> c)
 *
 * To determine which value maps back into 0, we compute:
 *
 *    (-c*minv(a,m)) % m
 *
 * and thus we find that the congruence generator would also normally map:
 *
 *    10239951819489363767 ==> 0
 *
 * To overcome this, and preserve the 1-to-1 and onto map, we force:
 *
 *    0 ==> 0
 *    10239951819489363767 ==> 1363042948800878693
 *
 * To repeat, this function converts a values into a seed value.  With the
 * except of 'seed == 0', every value is mapped into a unique seed value.
 * This mapping need not be complex, random or secure.  All we attempt
 * to do here is to allow humans who pick small or successive seed values
 * to obtain reasonably different sequences from the generators below.
 *
 * NOTE: This is NOT a random number generator.  This function is intended
 *     to be used internally by sa55rand() and sshufrand().
 */
define
randreseed64(seed)
{
    local ret;            /* return value */
    local a;            /* generator 0 multiplier */
    local c;            /* generator 0 adder */

    /* firewall */
    if (!isint(seed)) {
    quit "bad args: seed must be an integer";
    }
    if (seed < 0) {
    quit "bad arg: seed < 0 is reserved for future use";
    }

    /* if seed is 0, we will return 0 */
    if (seed == 0) {
    return (0);
    }

    /*
     * process the low 64 bits of the seed
     */
    a = 6316878969928993981;
    c = 1363042948800878693;
    ret = (((seed & cry_mask)*a + c) & cry_mask);

    /*
     * Normally, the above equation would map:
     *
     *        f(0) == 1363042948800878693
     *        f(10239951819489363767) == 0
     *
     * However, we have already forced f(0) == 0.  To preserve the
     * 1-to-1 and onto map property, we force:
     *
     *        f(10239951819489363767) ==> 1363042948800878693
     */
    if (ret == 0) {
    ret = c;
    }

    /* return the scrambled value */
    return (ret);
}


/*
 * Initial read execution code
 */
cry_seed = srand(0);        /* pre-initialize the tables */


global lib_debug;
if (lib_debug >= 0) {
    print "ver: 10.5 06:19:34 5/9/93";
    print "shufrand() defined";
    print "sshufrand(seed) defined";
    print "rand([a, [b]]) defined";
    print "srand(seed) defined";
    print "cryrand([a, [b]]) defined";
    print "scryrand([seed, [len1, len2]]) defined";
    print "scryrand(seed, ip, iq, ir) defined";
    print "random([a, [b]]) defined";
    print "srandom(seed) defined";
    print "obj cryobj defined";
    print "randstate([cryobj | 0]) defined";
    print "nxtprime(n, [val, modulus]) defined";
}